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We investigate the characteristic properties of self-sustained MRI turbulence in 
low-ionized proto-planetary disks. We study the transition regime between active and 
dead-zone, performing 3D global non-ideal MHD simulations of stratified disk cov- 
ering range of magnetic Reynolds number between 2700 < R m < 6600. We found 
converged and saturated MRI turbulence for R m ^ 5000 with a strength of a ss ~ 0.01. 
Below R m < 5000 the MRI starts to decay at the midplane, having Elsasser numbers 
below one. We find a transition regime between 3300 % R m ^ 5000 where the MRI 
turbulence is still sustained but damped. At around R m ^ 3000 the MRI turbulence 
decays but could reestablished due to the accumulation of toroidal magnetic field or 
the radial transport of magnetic field from the active region. Below R m < 3000 the 
MRI cannot be sustained and is decaying. Here hydro-dynamical motions, like den- 
sity waves dominate. We observe long-living anti-cyclonic vortices in the transition 
between dead-zone and active zone. 



Subject headings: accretion discs, magnetohydrodynamics (MHD), dynamo 



Introduction 



The magneto-rotational instability (MRI) i s a candidate for driving turbu l ence and enabling 



the ac cretion of matter onto the central object (|Balbus & Hawleylll99il 1 19981 : lHawley & Balbus 



199 1|) . Dependent on the ionization degree, the MRI generated turbule nce will b e reduced, down t o 



a low turbulence regime, called the dead zone (IBlaes & Balbuslll994l ; |jinlll996l : ISano et al. 



Various studies showed tha t a certain level of resistivity suppresses MRI act i vity (Fleming et al. 



2000; 
2010l : 



Turner & Sano 



Sano & Stone 2002aJJbk Fleming & Stone! 120031 : llnutsuka & Sand 120051 : iTurner et al. 



2000). 



2007, 



2008J). One of the most important dimensionless number, characterizing the 



coupling between gas and magnetic fields, is the magnetic Reynolds number R m which relates the 
timescale of magnetic advection to magnetic diffusion. We consider here the Ohmic diffusion term 
which is most important at the midplane. Up to now, there is no clear pres cription for which val- 
ues of R m the MRI turbulence is sustained in disks. Fleming et al.l (|2000|) found sustained MRI 
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with a zero-net flux magnetic field for R m > 10 4 . A recent study by ISimon et al.1 (1201 using 
stratified local box simulations, investigate the turbulence level for different values of the magnetic 
Reynolds and Prandtl number. They found a so-called low state, a state where turbulence is partly 
suppressed, but sustained due to a dynamo process. In addition, they predict a critical Reynolds 
numbers RZ lt in the rang e between 3200 < R% lt < 6000. A similar region was investigated by 
Qishi & Mac Lowl (1201 II) in which they found R™' ~ 3000. This critical Reynolds number is im- 



po rtant to mode l the su rface density of active layers in proto-planetary disk as it was recently done 
bv lMartin etaD d2012h . 

In our study, we will search for the critical magnetic Reynolds number in global zero-net flux 
stratified MRI simulations. Here the MRI turbulence criterion, Elsasser number v\ /(r/Q.), should 
become unity. In contrast to the magnetic Reynolds number, the Elsasser number gives one clear 
threshold independent of the magnetic geometry or the stratification. We w ill also investigate 
the hy drodynamical motions which become important in the dead-zone region (lOishi & Mac Low 



20091) . We concentrate on the magnetic Reynolds number regime proposed bv lSimon et al.l (1201 II) 



and 



Qishi & Mac Lowl(|201 lb . For our simulations we use only explicit resistivity. 



Qishi & Mac Low 



(1201 lb found out that well ionized MRI turbulence scales independently of Pm if Rm > Rm cr " . In 
addition, as the molecular viscosity is very small in proto-planetary disks we expect Prandtl num- 
bers of Pm « 1 and we focus on this low Prandtl number regime Q. In this paper we will first 
describe the numerical and physical setup. Then we will present the results, discussion and the 
conclusion. 



2. Setup 

The initial conditions for density, pressure and azimuthal velocity follow hydrostatic equilib- 
rium. We set the density p to 

„ v? /sin(0) - 1\ 

with po = 1.0, the scale height to radius H/R = c = 0.07, R = rsin(#). The pressure follows 
locally an isothermal equation of state: P = c]p with the sound speed c s = Co/ The azimuthal 
velocity is set to 

r \ sin(60 

The initial velocities V r and Vq are set to a white noise perturbation amplitude of V'" g " = 10 _4 Cj. 
We start the simulation with a pure toroidal magnetic seed field with constant plasma beta {3 = 



V*= ■ ( 2 ) 



Here the hydrodynamic viscosity is small compared to the magnetic diffusivity v «r) 
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2P/B 2 = 25. To obtain a range of magnetic Reynolds number we keep the magnetic dissipation 
value constant in the disk. We use three different values of 77, 771 = 2 • 10~ 6 AU 2 /yr, r/ 2 = 2.6 • 
\Q- 6 AU 2 /yr and 773 = 3.2 • \0' 6 AU 2 /yr. 



CjH 



= 2450 



\0.0l) \\AU! '{ 



2 ■ 10- 6 



(3) 



To estimate the numerical magnetic Reynolds number having MRI turbulence simulations we run 
also ideal MHD simulations with different resolutions as reference. The radial domain extends 
from 1 to 10 AU. The 9 domain covers + 4.3 disk scale heights, or 9 = n/2 + 0.3. For the azimuthal 
domain we use 2n for the L models and n/2 for the H models. We use a uniform grid in spherical 
coordinates. Models L have a resolution of N r = 384, N g = 192 , = 768 and 7Y = 384 for the 
H models. 

All models resolve the radial scale height with 9 to 22 grid cells for the inner and outer radius. 
The vertical scale height is resolved by 22 grid points. In models L the azimuthal scale height is 
resolved by 9 grid cells. The H models have a higher resolution of 17 per scale height in the az- 
imuth. Their are calcu lated with the FARGO MHD to reduce even more the numerical dissipation 
(|Mignone et al.ll2012|) . The simulation models are summarized in table 1. We note that model L l 
is special. Here the numerical dissipation cannot be neglected. By comparing with the results of H 
models, the L 1 model shows magnetic Reynolds number below the value used in H 3 . This model 
establishes a large dead-zone region. Here hydro-dynamical motions become important. Buffer 
zones extent from 1 to 2 AU as well as from 9 to 10 AU. In the buffer zones we use a linearly 
increasing resistivity (up to 77 = 10~ 3 ) reaching the boundary. This damps the magnetic field fluc- 
tuations and suppresses boundary interactions. For our analysis we use the range between 3 and 8 
AU, which is not affected by the buffer zones. Our outflow boundary condition projects the radial 
gradients in density, pressure and azimuthal velocity into the radial boundary and the vertical gra- 
dients in density and pressure at the 9 boundary. For all runs we employ the se cond order scheme 
in the PLUTO code with the HLLD Riemann solver (IMiyoshi & Kusanoll2005l) . piece-wise linear 
reconstruction and 2 nd order Runge Kutta time integration. We treat the induction equation with 
th e "Constrained Transport " (CT) method in combination with the upwind CT method described 



in 



setup can be found in 



Gardiner & Stond ( 2005b. using expli cit resistivity. A more detailed description of the physical 



Flock etal.1 (1201 11) 



3. Results 

According to Eq. 3, we obtain a specific value of magnetic Reynolds number at each radius, 
resulting in a specific value of turbulence. To compare the results at different specific radii the 
timescale of radial mixing is important. Due to the low magnetic dissipation, the timescale of 



-4- 



0.020 r 



OL 



Ideal FARGO 



- " y Ideal 
. ~* Tjldeal FARGO 

X tl 

0.015 - nHl 
- H 



Ideal MHD 

1: 

[] _ 



0.010 



0.005 - 



o.ooo L 



3000 



4000 



5000 



6000 



7000 



o.i8 r 



0.16 - 



. t Ideal 



0.141 

: - h 

0.12 - 



0.10- 
0.08 - 
0.06 - 

0.04 [ 

3000 



Ideal FARGO 
Ideal FARGO 



4000 



5000 



6000 



7000 



0.006 r 



0.004 



a* 



0.002 



Ideal FARGO 



_ t Ideal 

A L 

Tjldeal FARGO 

□ H 
H 



0.000 [ 



3000 



4000 



5000 



6000 



7000 



Fig. 1 . — Top: a S s value plotted for different magnetic Reynolds number R m for a time average of 
25 to 60 local orbits. Middle: Turbulent RMS velocity in units of the local sound speed. Bottom: 
Dynamo of the northern hemisphere. The L models are plotted in blue the H models are 
plotted in red. The error bars shows the standard deviation in the time average. Values for the 
corresponding ideal MHD runs are plotted at Rm = 6900 as reference. The accretion stress and the 
turbulent velocities saturate at magnetic Reynolds numbers above 5000. 
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magnetic diffusion is very long compared to the turbulent mixing timescale. The maximum radial 
mixing scale for an alpha value of 0.01 after 1000 years at the outer radius at 8 AU is around 
5x = y/2a ss c s Ht ~ 0.5AU. For our analysis we use specific radial positions which are separated 
enough. In chapter 3.5 we will investigate how the radial transport of magnetic field effects the 
local evolution. The a ss value at a specific radial position is calculated with 



ass(r) i — j& — >■ (4) 

The integral is done for each radius separately J r dV = J^ r J e f r 2 sin 6drd6d(p with r t and r r is 
the left and right radial cell boundary. The a ss is mass weighted to recover the correct total value. 
In chapter 3.1 we compare results having the same local time period. Here we average a ss (r), 
the turbulent velocities Vrms - ^{8v r ) 2 + {5v e ) 2 + {8v<p) 2 and the dynamo between 25 and 60 
local orbits. This ensures that we compare the same dynamical evolution at a specific radius and 
Reynolds number. The space average of Vrms and is done at the midplane region between 
0-1.5 scale heights. In chapter 3.2 we concentrate on the longterm evolution for different values 
of R m . There we include the two MRI criteria, the Elsasser number 

prjQ 

as well as the factor Q 




\5pc 2 RA(P' { ' 

The Q factor is the ratio of the MRI fas test-growing azim uthal wavelength to the azimuthal cell 



size. The value should be larger than 8 (|Flock et al.ll2010h . In chapter 3.3 we concentrate on the 



run L 1 having the largest dead-zone. Here we present results of the dominant hydro-dynamical 
motions. 



3.1. Time averaged statistics 

Due to the change of the rotation period Q with radius, the comparison of the turbulence 
statistics at different radii is limited to the number of rotations at the outer radius. We therefore 
present the comparison for the different Reynolds numbers between a time average of 25 - 60 
local orbits. We note that in this period the initial net magnetic flux has already vanished and 
we have a zero-net flux MRI turbulence. We note also that the results from model L 1 should 
be excluded. They indicate that here the explicit dissipation is not resolved and the numerical 
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dissipation dominates. We still overplot the results obtained by this lower resolution (L models, 
Fig. 1 blue color) to estimate the total magnetic Reynolds number. This model will be still useful 
as it shows very similar conditions as present in dead-zones. 

Fig. 1 combines the results of a ss , the turbulent RMS velocity V RM s and the dynamo value 
as a function of the magnetic Reynolds number. The results show for all turbulent quantities a 
saturation for magnetic Reynolds numbers around 5000 and above. Below R m = 5000 the MRI 
is damped and the turbulence decreases. We observe a slightly higher turbulence level for high 
magnetic Reynolds numbers compared to fully ideal MHD run. 

There is a saturation of ass at around 0.015 value for magnetic Reynolds numbers greater than 
5000, see Fig. 1 top. Below R m < 5000, a S s drops down to 0.005 at around R m ~ 3300. Even one 
could see a flattening at around R m = 3300 the values are still dominated by Maxwell stress which 
indicates that the MRI is still operating. The turbulent velocity scales roug hly with the square 
root o f V R ms ~ V^ c f found in recent local box simulations of dead-zones (|Okuzumi & Hirose 



201 lb . We observe a saturation around V RM s = 0.1 5c s (Fig. 1, middle). Below R m < 5000, 



V R ms drops down to 0.09 at around R m ~ 3000. We note again that these values represent the 
turbulent motions at the midplane, in the corona the turbulent velocity increases also in the dead- 
zone models. The dynamo action could be expressed by showing the correlatio ns EMF'^ = cumBj, 
with EMF'^ = v' r B' e - v' e B' r . We follow the analysis described in chapter 3.7 in (|Flock et al.ll2012h 
and a^ is normalized over c s . The dynamo action in the northern hemisphere is shown in Fig. 1, 
bottom, a^ saturates at around 0.004 and decreases only down to 0.002 at R m = 3000. In all runs 
the sign of in the southern hemisphere is negative with similar amplitude as the corresponding 
value in the northern hemisphere. 

The results of the ideal MHD runs show convergence. Here the FARGO MHD method plays 
an important role in decreasing numerical dissipation as it is presented in Mignone et al. 2012, 
A&A accepted. a S s converge around 0.01. We note again that the results presented in Fig. 1 
reflect mainly the turbulence level at the midplane. In addition we see here no decay of the MRI 
turbulence due to the short time period. In the next chapter we will concentrate on the longterm 
evolution over height for specific Reynolds number. 



3.2. Longterm evolution 

The previous chapter has shown the transition to saturated turbulence in well ionized disk 
regions. In this chapter we focus on the longterm evolution, including the turbulence over height. 
For each model we perform the analysis at a specific Reynolds number which are located in the 
middle of the domain. Here, H l has R m = 5500, R m = 4300 for model H 2 and R m = 3300, 3000 for 
model H 3 . In model L 1 we expect to have a Reynolds number around R m < 3000 by comparison 
with model H 3 . 
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We plot the time evolution of (255 in Fig. 2. The analysis show that both ideal runs as well as the 
resistive runs down to R m = 4300 (H l and H 2 ), show a steady turbulent evolution with a ss ~ 0.01 . 
The models H 3 and L 1 with R m < 3000 show a steady decrease of the turbulence until they oscillate 
in the range between 10~ 4 - 10 3 . Due to the decrease we note that the results of the models H 3 and 
L 1 plotted in Fig. 1 depend on the averaged time period. For these models, an extended simulation 
runtime was needed. We mark the regions with dominating Reynolds stress in thick (blue and 
green solid lines). The plots show only a short-term domination which indicate that the total 
integral over height suggest still a operating MHD turbulence. The configuration of the Maxwell 
stress cxm = B' R B'JAnpc 2 s over height is presented in Fig. 3. After around 200 inner orbits (years) 



the initial net azimuthal flux is lost and the fl ux start to oscillate around zero (IFromang & Nelson 



2006l : lBeckwith et aU201 ll : Flock et aU201 ll) . For a magnetic Reynolds number above R m = 5500 



the ionization is well enough to sustain the MRI turbulence at each height. In the midplane region 
we found a saturation of a ss ~ 0.01 at around R m > 5000. In model H 1 , Fig. 3, top, the Elsasser 
number is between 1 and 10 in the midplane region. If the ionization is decreased below R m = 5000 
the Elsasser number reaches unity at the midplane (see also Fig. 4). In model H 2 , Fig. 3, second 
from top, there is a small region at the midplane with 0. 1 < A < 1 . The total ass is still around 0.01 
(compare Fig. 2) but with higher fluctuations. Locally the magnetic turbulence starts to vanish at 
the midplane. At around R m = 3300 (model H 3 ) the total a M gets affected and decreases (compare 
Fig. 2). In model H 3 the Elsasser number A = 1 point is around 1 scale height. At the midplane 
A drops below 0.1. We confi rmed the A = 0.1 region equivalent with the X res = H definition 
by lOkuzumi & Hirosd (|201lb . Fig. 1 therein. At this point also the azimuthal MRI wavelength 
becomes unresolved (Q < 8, red solid line). At the midplane the Reynolds stress dominate over 
the Maxwell stress (Fig. 3, yellow solid line). At R m ~ 3000, the MRI gets strongly damped. For 
model L 1 and model H 3 after around 500 years, the poloidal magnetic fields decreases until the 
Elsasser number drops below 0.1 as well as the azimuthal MRI wavelength becomes unresolved 
(Fig. 3, red solid line). If the turbulence is able to survive depends also on the surrounding layers. 
Even we choose a constant resistivity, the turbulence in the runs H 3 and L 1 is similar as they are 



present in dead- zones. Turbulent u pper layers of the disk are usually active (|Turner & Sanoll2008 



Dzyurkevich et al .1120 10c LB ail 1201 ll) . In this active layers the turbulent mixing is also stronger and 



active channels can reach into the dead-zone. Fig. 4 presents a 2D r—6 slice of radial magnetic field 
overplotted with the A = 1 line for model H l and H 2 after 1000 years. At around R m = 5000 the 
dissipation starts to damp the MRI at the midplane while there is still MRI turbulence in the corona 
of the disk. A similar shape could be expected at the outer edge of dead-zones (Dzyurkevich et al. 
2012, ApJ subm.). 

We summarize that for magnetic Reynolds numbers above R m > 5000 we observe a saturated and 
converged MRI turbulence in zero-net flux stratified simulations down to the midplane region. The 
Maxwell stress dominates in each region the Reynolds stress. Below R m < 5000 the turbulence 
starts to decay at the midplane. The Elsasser number drops below 1 and there regions with low 
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magnetic fields and dominating Reynolds stress. At around R m = 3000, the Elsasser number 
drops below 0. 1 at the midplane. Long-period oscillations becomes visible with MRI activation 
and decays. Below R m < 3000 the MRI cannot operate. Here, Reynolds stress dominates over 
the Maxwell stress. Dependent on the surroundings, the upper layers are still active having MRI 
turbulence. The results show that for magnetic Reynolds numbers down to R m = 4300, the MRI 
driven by a zero-net flux azimuthal magnetic field, can sustained the turbulence with ass - 0.01. 
Models H 3 and L 1 with R m < 3000 show a decrease of the turbulence and eventually a dominating 
Reynolds stress. Fro m this results we conclude that the critical magnetic Reynolds number should 
around R%" < 3000. lOishi & Mac Lowl (|201ll) also found a critical magnetic Reynolds number 
around R m ~ 3000 in local box simulations. In the next chapter we concentrate on the run L 1 which 
present similar conditions as they are present in dead-zones. 



3.3. Dead-zone oscillation 

Before we move to the hydro-dynamical motions in dead-zones we want to investigate the 
long-term oscillations of the Maxwell stress we observ ed in the previous c hapter. This oscillations 



could be connected to some sort of dynamo process. iSimon et al.1 (|2011|) presented in local box 
simulations a sudden transition from a 'low' turbulent state to a 'high' turbulent state. We observe 
in model H 3 and L 1 for R m ~ 3000 a similar effect. The results are combined in Fig 5 showing 
the mean toroidal magnetic field overplotted with the turbulent EMF'^ (times 1000), the Elsasser 
number and the Q value. For this analysis we choose the southern hemisphere between and 1 .5 
scale heights. After around 200 years the initial mean toroidal magnetic field (black solid line) 
vanishes and starts to oscillate. There is a clear correlation between the turbulent EMF (red dotted 
line) and the mean field which indicates a working aO. dynamo. In both models, after around 500 
years the MRI from azimuthal magnetic fields becomes unresolved (Q < 8 blue dotted line below 
black dotted line). Here the Elsasser number drops below unity which still allows for damped 
MRI growth. After 800 years the Elsasser number drops below 0.1. At this stage there is no MRI 
working and we are in the dead-zone stage. Here also the oscillations of mean field and turbulent 
EMF stop. In model H 3 the MRI switch on again after less then 200 years (27 local orbits) with 
a larger oscillation period. In model L 1 we observe a late switch on at 1100 years for only 150 
years (20 local orbits). Afterwards the Elsasser number reaches nearly 10~ 3 and we don't expect a 
relaunch again. The sudden increase of A is due to a sudden increase of poloidal magnetic field. 
This sudden increase could be connected with the accumulation of toroidal magnetic fields (black 
solid line in Fig. 5, bottom) in the dead-zone as well as due to radial transport of magnetic field 
(see chapter 3.5). In Fig. 6, we present an azimuthal slice of toroidal magnetic field in the r — 6 
plane after 1500 years for model L 1 . In the dead-zone (A < 0.1) we observe locally a Q factor of 8 
and above which correspond to a plasma beta below 50, still the dissipation is too high to relaunch 
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Model name 


N,-Nf>N* 

* * T * u * (p 


Ar : A6 : Aphi 


FARGO 


Rm 


j^Ideal 


384x192x768 


1-10: 0.6 


In 


NO 


Ideal 


jJdealFARGO 


384x192x768 


1-10: 0.6 


In 


YES 


Ideal 


j_jIdealFARGO 


384x192x384 


1-10: 0.6 


n/1 


YES 


Ideal 


H 1 


384x192x384 


1-10: 0.6 


n/2 


YES 


4300-6300 


H 2 


384x192x384 


1-10: 0.6 


n/1 


YES 


3300-4800 


H 3 


384x192x384 


1-10: 0.6 


n/1 


YES 


2700-4300 


L 1 


384x192x768 


1-10: 0.6 


In 


NO 


<3000* 



Table 1 : From left to right: model name, resolution, domain size, FARGO-MHD, range of magnetic 
Reynolds number. In model L l the used explicit resistivity is not resolved and we estimated the 
magnetic Reynolds number. 
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Fig. 2. — a ss time evolution for different models. The value is taken in the middle of the domains 
for H l at R m = 5500, for model H 2 at R m = 4300 and R m = 3300 for model H 3 . For model L 1 we 
assume R m < 3000. The ideal MHD models as well as the resistive models down to R m = 4300 
show a values around 0.01. We mark the time with dominating Reynolds stress with a thick line. 
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Fig. 3. — Evolution of the Maxwell stress % over height. For model H l at R m = 5500 (top), H 2 
at R m = 4300 (second from top), H 3 at R m = 3300 (third from top) and L 1 at R m < 3000 (bottom). 
The black solid lines show the Elsasser numbers 10,1, and 0.1. The red solid line shows 2 = 8. 
The vellnw line shows the re.pinn with Hnminatinp Revnnlds stress 
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Fig. 4. — Contour plot of radial magnetic field in the r - 8 plane after 1000 years overplotted with 
the Elsasser number A = 1 line. In both models, H x (top) and H 2 (bottom) the critical value of R m 
of resolved MRI at the midplane region is around R'j"' = 5000. 




Fig. 5. — Time evolution of mean azimuthal field (black solid line), turbulent EMF'^ (red dotted 
line), A (blue solid line) and the Q value (blue dotted line) for model H 3 at R m = 3000 (3.8 AU) 
(top), and for model L 1 at 3.8 AU (bottom). Values are calculated in the southern hemisphere 
between (0 - 1.5 scale heights). The black dotted line shows the value of Q = 8. In model L 1 
we observe a relaunch of MRI between 1200 and 1300 years in the dead-zone. For model H 3 we 
observe only short-period switch off between 800 and 900 years. 
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the MRI growth by azimuthal fields. In general the line of resolved azimuthal field Q < 8 matches 
the one of decaying MRI A < 0.1 (compare Fig. 3). 



3.4. Hydrodynamical motions 

This dead-zone region is also interesting in terms of the hydro-dynamical motions. A close 
looks at the radial velocity in the r - <p midplane in Fig. 7 reveals velocity amplitudes around 
0.2c s . On the first look these hydro-dynamical waves appear to be linear waves with a single 
mode. But a closer look at the Fourier spectra of the radial velocity reveals a more detailed picture. 
We calculated the spectra in the dead-zone region between 4 and 5 AU, time averaged between 
1300 and 1500 years in model L 1 . In Fig. 8, we compare this spectra with the one obtained 
in fully ionized disk from model L Ideal FARG °. In the dead-zone there is a peak at m = 6. The 
turbulence at lower scales (green solid line) is around 1 order of magnitude below the value in 
the fully MRI turbulent region. Beside the density waves we observe anti-cyclonic vortices (Fig. 

7 between 6 and 7 AU) which create large extended spiral arms. They are produced at around 

8 AU where we have a positive density slope due to the buffer zones. We calculated the relative 
vorticity to = ((V x V)g - (V X V Kep i er ) e )/(V x v Kep ier)e between -0.8 and -0.5 in the large vortex. 
The growth of vortices at the border of the dead zone can be due to the Rossby wave instabi lity as 



proposed by lVarniere & Tagged (I2006|) and recently investigated by lLyra & Mac Low! (120 121) . The 
vortex structure is presented in Fig. 9. The vertical extension is around +2 SH, the radial extension 
is 2H r ~ IAU and the azimuthal one is around 10//^ ~ AAU . The vortex is dragged by the 
hydro-dynamical surroundings which gives him the concave (Fig. 9 , top, dominant radial inward 
velocity), convex (Fig. 9, bottom, V'^ > ) shape respectively. The magnetic fields (Fig. 9, 
bottom black vectors) are not present inside the vortex but in the surrounding layers. We measure 
a plasma beta of ~ 10 5 inside the vortex compared to ~ 100 - 1000 in the surrounding region. 



3.5. Radial mixing 

We already mentioned the effect of radial transport due to turbulent mixing. In this subsection 
we want to focus on the importance of radial transport of magnetic fields. We calculate the diver- 
gence of the radial Poynting flux J V • S r dt integrated over time for the total 6 and <p domain at 
A.6AU using one scale height in radius. We integrate the divergence of the radial Poynting flux over 
every output (dt = 0.1 local orbits) and compare it with the total magnetic energy. Fig. 10 shows 
how much magnetic energy is transported at each output (0.1 local orbits). In the fully turbulent 
layers the net radial transport is around zero with fluctuations of around 5 %. In the transition zone 
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Fig. 6. — Contour plot of azimuthal magnetic field in the r - 9 plane after 1500 years overplotted 
with the Elsasser number A = 0.1 red solid line and the 2 = 8 green solid line. We observe in the 
dead-zone relative strong (resolved) toroidal magnetic fields. 
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Fig. 7. — Contour plot of radial velocity in the r-(f> midplane after 1500 years in units of the sound 
speed. The dominating Reynolds stress in the dead-zone is produced by linear waves. 
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Fig. 8. — Fourier spectra of radial velocity over azimuthal wavenumber at the midplane compared 
for the active zone (model Jj deal FARG ° anc [ the dead-zone (model L 1 ). The density wave peak 
appears in both models at m = 6. We normalize with the amplitude at Vj deal (6) = 0.25c s . 
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Fig. 9. — Top: Contour plot inside the vortex of turbulent azimuthal velocity in the r — 6 plane 
overplotted with the velocity vectors. Bottom: Contour plot inside the vortex of turbulent radial 
velocity in the (p-6 plane overplotted with the magnetic field vectors. Color table present the units 
of the sound speed. 
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and the dead-zone the radial transport becomes more important and these fluctuations reach peaks 
of 10 to 20%. 



4. Discussion 



4.1. Turbulence between active and dead-zones 



In our work we showed that the MRI turbulence saturates around R m = 5000 with a ss ~ 0.01 . 
Here the Elsasser number is around A > 10. This turbulence level is sustained for Elsasser number 



values of 1 < A < 10. Stratified converged ideal MHD local simulations (|Davis et al.l l2010h 
presented similar values of ass- Still our results indicate a higher turbulence level using small 
explicit resistivity. High resolution simulations using explicit resistivity are needed to confirm 
the convergence of the accretion st ress. In addition, there are mean-field mechanism, which can 
temporal increase the value of a S s (|Flock et al.ll2012l) . 

For magnetic Reynolds numbers around the critical value R m ~ 3000 we observe long-period 
oscillations of the accretion stres s as well as irregul ar switch off and on of MRI activity. Such 
oscillations were also reported by ISimon et al.l (|201 lb in local box simulations. We see indications 
that these oscillations are trigger by some dynamo process as locally there are accumulations of 
strong mean toroidal fields in the dead-zone. Also radial transport of magnetic field could here play 
a role. In such configurations the corona is still active while in the midplane the MRI is switched 
of. Here the Reynolds stress dominate over the Maxwell stress at A < 0.1. A similar classification 
of turbulent regions in pro to-planetary disks was done by comparing different heights, see also in 
Okuzumi & Hirosel (boi lb . Fig. 1. 



4.2. Numerical, explicit and turbulent dissipation 

In such turbulent simulations there are 3 kinds of important dissipation sources. The numer- 
ical dissipation, due to the finite grid size, the explicit resistivity and the turbulent dissipation. In 
model L 1 the numerical dissipation plays an important role. In addition one has to note that using a 
uniform grid, the numerical dissipation will be larger at the inner part then in the outer part. Here a 
logarithmic increasing grid would help. Then there is of course the explicit dissipation which is in- 
cluded in the induction equation in the code. Here one should be sure that the explicit one is higher 
than the numerical one. The turbulent dissipation scales with the level of t he turbulence. Here it is 



still un clear how this dissipation process works in detail, see also Fig. 9 in lFromang & Papaloizou 



(2007) 
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Fig. 10. — Time integral of divergence of radial Poynting flux over total magnetic energy for model 
H Y (R m = 5500, blue solid line), H 2 (R m = 4300, red dotted line), H 3 (R m = 3300, green dashed 
line) and L 1 (R m < 3000, black dash dotted line). In the transition region the radial mixing becomes 
more important. 
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5. Summary 

We performed 3D global ideal and non-ideal MHD simulations to study MRI turbulence in 
low-ionized proto-planetary disks using an initial toroidal magnetic field for a range of magnetic 
Reynolds numbers. With our global simulations we are able to investigate the transition regime 
which is present between the active and the dead-zone in proto-planetary disks. We define 3 dif- 
ferent disk regimes dependent on the magnetic Reynolds numbers: 

• Above a Reynolds number of R m <: 5000. Here MRI turbulence saturates and is sustained. 
For this region we find steady and converged ass values around 0.01 in ideal and non-ideal 
simulations. The Maxwell stress dominates down to the midplane region. The turbulent 
velocities reach maximum values around V RM s = 0.1 5c s . The aQ, dynamo is operating and 
we find = 0.004 with a positive sign in the northern hemisphere. The Elsasser number 
A stays above 1 in the midplane. 

• Reynolds number between R m = 3000 - 5000. Here, the MRI starts to switch of at the 
midplane region. The Elsasser number A drops below 1. Simulations with 3300 < R m < 
5000 show a still sustained turbulence, supported by radial transport of magnetic fields. For 
R m = 3000 the MRI turbulence starts to switch off. We expect the critical magnetic Reynolds 
number around R"j" ^ 3000 and below. 

• Reynolds number around R m % 3000 and below. Here, there is no MRI at the midplane 
anymore. The Elsasser number A drops below 0.1. The turbulence is dominated by the 
Reynolds stress. We observe long-period oscillations of MRI activity and MRI decay. There 
is an accumulation of toroidal magnetic fields below /3 < 50 in the dead-zone. The turbulence 
at the midplane is supported by active channels in the corona which pump sound wave into 
the midplane region. The velocity spectra at the midplane reveals a drop of one order of 
magnitude at the scale of H in the dead-zone compared to the fully ionized turbulent region. 
We observe long lived anti-cyclonic vortices in the transition regime, creating spiral arms in 
the dead-zone. Magnetic fields are not present in the inner part of the vortex. 
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Natalia Dzyurkevich for her suggestions and comments during this work. We thank also Neal 
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from the European Research Council under the European Union's Seventh Framework Programme 
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